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Abstract 

We present a new method for the numerical evaluation of arbitrary loop integrals in dimensional 
regularization. We first derive Mellin-Barnes integral representations and apply an algorithmic 
technique, based on the Cauchy theorem, to extract the divergent parts in the e — > limit. We 
then perform an e-expansion and evaluate the integral coefficients of the expansion numerically. The 
method yields stable results in physical kinematic regions avoiding intricate analytic continuations. 
It can also be applied to evaluate both scalar and tensor integrals without employing reduction 
methods. We demonstrate our method with specific examples of infrared divergent integrals with 
many kinematic scales, such as two-loop and three-loop box integrals and tensor integrals of rank 
six for the one-loop hexagon topology. 
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I. INTRODUCTION 



Perturbative methods are indispensable in order to establish consistent theories of particle 
interactions, and to predict quantitatively their experimental manifestations. The anticipa- 
tion of new phenomena in modern experiments and theoretical extensions of the Standard 
Model, requires cross-sections for complicated processes. This has driven a remarkable 
progress in the development of new computational methods. At the one-loop level, the 
calculation of cross-sections with five external legs is gradually becoming a routine activity 
(e.g. (y). At two- loops, there have been recent successful computations of amplitudes with 
four external legs and up to three kinematic scales (e.g. 0). At three-loops and beyond, 
amplitudes with up to one parametric variable have also been computed (e.g. 0]). 

Many new processes with higher final-state multiplicity, number of loops, and kinematic 
scales, have been identified to be important at the TeV energy frontier. The aim of our 
paper is to provide a new method which can be used to compute loop amplitudes for such, 
more complicated, processes. 

Loop integrations are cumbersome due to the presence of infrared singularities. Loop 
integrals with many kinematic scales have, in addition, a complicated analytic structure 
with respect to their kinematic parameters. The tensor structure in gauge theories is also an 
issue, since it proliferates the number of terms. A general method for computing arbitrary 
loop integrals should extract their infrared and ultraviolet singularities, and treat simply 
kinematic discontinuities and threshold singularities. For practical applications, it should 
also be able to handle tensor integrals efficiently There is no method which addresses 
satisfactorily all these issues; known techniques can compute a limited number of amplitudes 
where some simplifications occur in special cases. 

Following a traditional approach to calculate loop amplitudes, one reduces the number 
of terms to a few master integrals and computes the latter analytically. One-loop integrals 
can be reduced using the classical method of Passarino and Veltman [4] . For generic multi- 
loop computations one can derive reduction identities from integration by parts [^] and the 
invariance of scalar integrals under Lorentz transformations jfj. There is a large variety of 
approaches for evaluating the master integrals analytically. For example, one can compute 
integrals with a simple singularity structure and a small number of kinematic scales by 
integrating directly their Feynman parameter representation. For more complicated cases 
one can use advanced techniques, such as the method of differential equations 0, 0, H| . 

This approach can fail, however, if we apply it to complicated processes. The reduction 
algebra is hard, and the expressions of amplitudes in terms of master integrals may have 
spurious singularities which hamper their numerical evaluation. The extraction of e poles 
in the master integrals, the evaluation of the coefficients of the e-expansion in terms of 
known analytic functions, and the analytic continuation of the latter to physically interesting 
kinematic regions are also involved. It is, thus, very well motivated to improve or replace 
the "traditional scheme" and to develop new automated methods. 

One-loop amplitudes can be entirely determined in four dimensions in terms of basic 
functions such as logarithms and polylogarithms that appear in the one-loop scalar box. New 
methods introduce sophisticated algorithms for the reduction to the basic functions; by either 
numerical or analytical techniques, they control the appearance of spurious singular terms 
and minimize the size of the intermediate expressions 0, HI 11 , 12^ l| 14]| • Recently, cross- 
sections for e + e~ — > 4 fermion processes at NLO were computed |l5| using such techniques. 
A method, which is inspired by techniques for real radiation at NLO, renders one- loop 
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graphs numerically integrable 0] with the introduction of universal subtraction terms. A 
different approach uses unitarity, dualities, and analyticity properties for the determination 
of one- loop amplitudes |l7| |. 

Beyond one-loop, there has been a significant progress in the automation of the reduction 
methods to master integrals [3, 1^, 2(|. The infrared structure of multi-loop integrals is sub- 



stantially more complex than at one-loop; the development of methods for their numerical 
evaluation is more difficult. Nevertheless, there is significant progress in this direction |2l| . 
A powerful numerical method for multi-loop calculations is the method of sector decompo- 
sition jiil . I23I . 24j : it simplifies recursively the singularities of Feynman parameterizations 



and allows a straightforward expansion in e. The method of sector decomposition has been 
very succesfully employed in the computation of several multi-leg integrals at one, two and 
three loops l22i 23] . This method has been introduced, recently, for the purely numer- 
ical evaluation of multi-loop amplitudes |26| . However, it is perplexing how to apply it for 
loop-integrals in non-Euclidean regions. 

In 1999, Smirnov j2?J and, soon later, Tausk j28[, introduced a new method for the 
evaluation of loop integrals. In their pioneering papers, Smirnov and Taus k 1271 l28j | computed 
analytically the first infrared divergent double box integrals. In [2^, l30j ] new two-loop 
integrals for 2 — > 2 massless processes were computed. In (3lJ the method was applied to 
double-box integrals with one additional mass-scale. The method was spectacularly applied 
in the computation of three-loop amplitudes [3^, HH in J\f = 4 supersymmetric Yang-Mills 
theory, shedding light to novel cross-order perturbative relations of the theory 33|, 34^ . 

The Smirnov- Tausk method is based on a few simple ideas. Starting from the Feynman 
parameterization of a loop integral we can derive a new representation in terms of a multiple 
complex contour integral. Such, Mellin-Barnes (MB), parameterizations have enabled com- 
picated loop calculations by using powerful methods for complex integration (36[. Smirnov 
and Tausk exploited a novel property of these representations. Infrared divergences localize 
on simple poles inside the complex integration volume. We can isolate the divergent pieces 
of the integral at e = 0, by using the Cauchy theorem. After subtracting the divergent 
residues, we can perform a Taylor expansion in e, and sum up the remaining infinite series 
of residues. Finally, we can work to derive analytic expressions for the infinite sums in the 
coefficients of the expansion in terms of logarithms, generalized polylogarithms |35j], and 
more complicated functions. 

The Smirnov- Tausk method is very powerful; however, it is laborious and intricate. The 
isolation of the divergent residues in multiple Mellin-Barnes integrals is convoluted. In ad- 
dition, it is difficult to identify infinite sums in terms of polylogarithm functions with known 
analyticity properties. As a consequence, the analytic continuation in physical kinematic 
regions is also involved. Due to these complications, the method has been applied to a few 
master integrals with a small number of kinematic scales. In this paper, we generalize the 
method to a broader spectrum of applications. 

As a first task, we automate the procedure for the isolation of the divergent residues 
at e — > 0. Smirnov j2?| and Tausk |28[ use different techniques for finding these residues. 
The approach of Smirnov is very intuitive, but daedal. We have found that the technique 
described by Tausk in Ref. 28[ resembles closely to a programmable algorithm. We have 
used it as a guide and we have written computer programs which subtract the 1/e poles in 
arbitrary Mellin-Barnes integrals. 

In our method, we avoid entirely the painstaking tasks of finding analytic expressions for 
infinite sums in terms of polylogarithms, and performing the analytic continuation in the 
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arguments of polylogarithms. Mellin-Barnes representations are valid in kinematic regions 
where loop integrals may be complex- valued. We have found that, in a broad spectrum 
of applications, it is simple to calculate the representations numerically. The only analytic 
continuation that is ever required is that of logarithms with a single kinematic scale as an 
argument. 

An important goal of our method is to calculate loop amplitudes in realistic gauge the- 
ories. We have found that tensor integrals and, furthermore, diagrams which belong to the 
same topology can be calculated collectively. As we will show, the integrand of a represen- 
tation for a scalar integral will be a product of Gamma functions and powers of kinematic 
invariants, while that of a generic tensor will be the same integrand as in the scalar integral 
multiplied by a polynomial in the integration variables. The evaluation of polynomials is 
fast in a numerical program; the computational cost for evaluating tensor integrals or loop 
diagrams is not significantly larger than evaluating scalar integrals. The only practical issue 
that we need to address, is the book-keeping of the various terms that contribute to the 
polynomial; we present an efficient solution of this problem here. 

In this paper, we apply our method to a number of examples. We, first, test our method 
in scalar and tensor integrals of the one-loop massless hexagon topology. The purpose of 
this computation is to introduce our method for tensor integrals and to demonstrate that we 
can tackle problems which are relevant in computations of physical amplitudes. We present 
here results for tensors through rank six, in both the Euclidean and the physical region for 
2 — > 4 processes. The numerical programs that we have constructed are suitable for the 
evaluation of the QCD amplitudes in four-jet production at hadron colliders. 

In the second set of examples, we compute scalar two and three-loop integrals which are 
known analytically in all kinematic regions: the massless planar Wh and cross 28] double- 
box, the massless double-box with one off-shell leg 31, 37, 38, Hjj, and the massless 
planar triple-box [3^ ■ They serve to cross-check our algorithms and to demonstrate that we 
can easily reproduce state-of-the-art computations. Our numerical results are in excellent 
agreement with the analytic expressions. We also present a number of new results that would 
require significant efforts for their computation with traditional approaches. We present, for 
the first time, double-box integrals with up to four kinematic parameters, and triple-box 
integrals with up to three kinematic parameters computed in all physical regions. 

In Section |nj we explain our technique for deriving Mellin-Barnes representations for 
loop integrals. In Section 11111 we describe our routines for performing an e-expansion of 
Mellin-Barnes representations for scalar integrals. We extend these results to loop integrals 
with tensor numerators in Section IIVI and, in Section we present our methods for the 
numerical integration. Section IVII is devoted to present our results for several integrals at 
the one, two and three loop level. Finally we present our conclusions. 



II. MELLIN BARNES REPRESENTATIONS 

We start with a brief discussion on the derivation of Mellin-Barnes representations for 
loop integrals from their Feynman parameterization. The construction of parameterizations 
is not unique, and various representations of the same integral may have quite different 
features. For example, they could have a different integral dimensionality, or they could 
be better or worse suited for numerical evaluation. A valid parameterization, however, can 
always be found. A pedagogical introduction to the topic of Mellin-Barnes representations 
is presented in Ref. j4l|. 
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As a concrete example we consider the one- loop box integral with two adjacent external 
legs off-shell and massless propagators, 

T 2m = f <Fk 1 / 1\ 

4 / . d A V\ A U2 A wz, A v& i \ > 

J %TX'2 /±i S\ 2 /I3 /I4 

with 

A 1 = k 2 + i0, 

A 2 = (k + Pl ) 2 + iO, 

A 3 = (k + Pl +p 2 ) 2 + iO, 

A 4 = {k + Pl + p 2 + p 3 ) 2 + i0, 

and v\=v\ = 0, (pi + P2) 2 = s, (p 2 + p 3 ) 2 = t, p\ = M 2 , (pi + p 2 + p 3 ) 2 = M 2 . The powers 
of the propagators and the dimension d = 4 — 2e are kept arbitrary and we will derive a 
Mellin-Barnes parameterization of the box-integral in this general case. In this section, we 
will seek values for {ui, d} where the representation is well defined and the integral is finite. 
In the next Section, we will explain the technique for the analytic continuation to the values 
of the parameters, such as{z/j = l,e = 0}, where the integral develops divergences. 
Our starting point is the Feynman parameterization of the box integral: 



j2m = ,_ V N 



r(;v-f) 



' 4 v T(y x )T(y 2 )T{v z )T{v±) 

r 1 dxidx^x^dx^Xi 1-1 x^ 2 ' 1 x 1 ^^ 1 x^ 1 S (xi + x 2 + x 3 + x A — 1) 

X / 9 o 7V—^ ' ^ ' 

J0 ( — sx\X 3 — tx 2 x 4 — Mfx^Xi — M 2 x^X\ — iO) 2 

with iV = v\ + v 2 + U3 + z/4. The main tool for obtaining the Mellin-Barnes representation 
of the above integral is the formula: 

1 1 rc+ioc r(—w)T(a + w) 

(A -\- A ) a = w/ ■ dwAfAr-" [ fvM ' (3) 
[A\ + A 2 ) ZmJc-too 1 [a) 

where the following conditions are satisfied: 

1. Ai t 2 are complex numbers with \args{Ai) — args (A 2 )\ < tt, 

2. the contour of integration is a straight line parallel to the imaginary axis separating 
the poles of T(—w) and T(a + w), i.e. — Re(a) < c < 0. 

We can use the identity of Eq. El to simplify the entangled denominator of Eq. |2l with the 
cost of introducing new integrations in the complex plane. Before we do so, we would like 
to point out an attractive feature of Eq. El The above two conditions guarantee that the 
integrand in the right hand side of Eq. El vanishes at infinity. We can then close the contour 
of integration at infinity and calculate the integral using the Cauchy theorem. If we close 
the contour to the side of the positive real-axis we obtain 



A^ T(a)n\ V A 2 



5 



This is the Taylor expansion of the left hand side of Eq. |3] in the region lAi/Aal < 1. 
Equivalently, by closing the contour of integration to the side of the negative real axis, we 
obtain a Taylor expansion in the complementary region IA2/-A1J < 1. This is a particularly 
useful property: representations of Feynman integrals which are derived by using the Mellin- 
Barnes decomposition are valid in all kinematic regions. In addition, A 12 can be complex; 
therefore, we can account for the infinitesimal imaginary part iO that is assigned to invariant 
masses and kinematic parameters. Eq. 3 may be recursively applied to denominators with 
more than two terms, yielding: 

1 rc+ioo 

\W m -l A—a—W\_ — ...—W m -l 



] rc+ioo 

-^n / d Wl ... dw m - X A^ . . . Al^A- 

(5) 



(A1+A2 + ... A m ) a (2m 

r(— wi) ■ ■ ■ r(-u? m _i)r(« + w 1 + . . . w. 



r (a) 



m) 



We are now ready to apply the decomposition of Eq. |5] to the denominator of Eq. El intro- 
ducing three contour integrations: 



Xl m (K d}) = (-1) N I f[ )S(x 1 +x 2 + x 3 + x 4 - 1) 

lr d 

x ^niy J dw i dw 2 dw 3F(-wi)~r(-w 2 )T(-w 3 )r(N - - + w 1 + w 2 + w 3 ) 

x (-x 3 x 4 M*) Wl (-x 4Xl Ml) W2 (-x 2 x 4 t) W3 (- XlX3 s)-2- N - w ^-^ , (6) 

where the contours must satisfy condition 2 above, i.e. they must be chosen in such a way 
that the arguments of the Gamma functions are all positive. It is now straightforward to 
integrate out the Feynman parameters, using 

To exchange the order of integrations, we must assume that the exponents of the Feynman 
parameters satisfy 

Re(oi) > 0. (8) 

With these conditions the Gamma functions on the right hand side of Eq[7]have arguments 
with positive real part as well. The final result is: 

(-l) N r 

ll m ({u t ,d}) = ^—i- J dw 1 dw 2 dw 3 r(-w 1 )V(-w 2 )V(-w 3 ) 

x f&l, Vz, V3, Vi, Wl, w 2) w 3 ) 

x (_ Ml 2 - i0) m (-M| - itf" 2 (-t - zOr (-s - zO)^^-" 123 , (9) 



where 



/(f 1, V2, V3, "A, Wl, W 2) W 3 ) = [Yl 



1 1 \ r (1/1234 - 5 + W 123 



[ T(ui) J T(d - 1/1234) 



xr (v 2 + w 3 ) r (1/4 + w 123 ) r I- - u 2U - w 13 \ r I- - v 124 - w 23 J (10) 
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and we have adopted the shorthand notation 



Vijk... = n + Vj + u k + . . . 

In general, Mellin-Barnes representations for a loop-integral are well defined for number 
of dimensions and powers of propagators which render the real parts of the arguments in 
all Gamma functions positive. For example, it is impossible to satisfy positivity in all 
Gamma function arguments for {ui = 1, d = 4}. This is in accordance with the expectation 
that the scalar box integral is infrared divergent in four dimensions. But, for example, the 
representation of Eq. [SJis well defined if we choose the contour Re(u>i) = Ref^) = Re(u> 3 ) = 
—0.2 and the parameters {z/j = l,d = 5.4}. In the next section we will detail the analytic 
continuation from values of parameters for which the Mellin-Barnes representation is well- 
defined to values for which the integral develops divergences. That is, a procedure which 
makes the divergences appear explicitly in the form of poles in the parameter space. 

From Eq. |§] we observe that it is simple to implement the analytic continuation of kine- 
matic parameters in a numerical evaluation, since they always appear in a factorized form, 
e.g. 

(-t-i0) w = e u,1 °g(-*- i0 ). 

We need the analytic continuation of simple logarithms; these are evaluated trivially in all 
kinematic regions: 

log(-*-i0)=log(|f|)-e(*)i7T. 

It is clear from the above example that we can derive a Mellin-Barnes representation for 
any integral starting from its Feynman parameterization. The form of the parameterization 
depends on the choice of Feynman parameters, the order in which the A{ terms appear in 
the right hand side of Eq. and the implementation of the delta function constraint for the 
Feynman parameters. This arbitrariness is more pronounced beyond one-loop. 

In Section HV1 we will derive representations for tensor one-loop integrals which simplify 
their numerical evaluation. We aim to find similarly simple representations for multi-loop 
tensor integrals as well. For this purpose, it is convenient to use the one-loop representations 
as building blocks, by employing a re-insertion technique [3(l|36j]. We explain this technique 
in our second example. 

We now derive a Mellin-Barnes representation for the massless planar double-box integral 
with two adjacent legs off-shell. The integral is defined: 

j2m = f d'k dH_ 1 (n) 

J Z7T2 Z7T2 AiA2A^A^A^AqA-j 



with 



A x = k 2 + i0, 

A 2 = (k + pif + iQ, 

A 3 = {k + Pl +p 2 ) 2 + i0, 

A 4 = (l + Pl +p 2 ) 2 + t0, 

A b = (I+P1+P2 + P3) 2 + 

A Q = l 2 + iO, 

A 7 = (k-l) 2 + iO, 
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and Pi = p\ = 0, (pi+p 2 ) 2 = s, (P2+P3) 2 =t,p% = Mf, (P1+P2+P3) 2 = Mf . For simplicity, 
we have set the powers of propagators to one. We perform, first, the fc-loop integral and 
derive the representation for this integral reading it off from the one-loop result of Eq. 

Jl m = I ^rx-hr / dw 1 dw 2 dw^{-w 1 )Y{-w 2 )Y{-w^) 

x/(l,l,l,l, Wl , U ; 2 , W 3)(-s)^ 4 - Wl23 (-^4) UJl (-A6r 2 (-^8)" 3 ; (12) 

the off-shell legs of the /c-loop box correspond to the A4, A$ propagators while the t invariant 
mass gives rise to a new propagator for the second integration: 

A 8 = (l + Pl ) 2 . (13) 

The propagators A4, A5, Aq, As form a new one- loop box with two adjacent legs off-shell, 
and we can read off, once again, the representation for this second integration using Eq. 
The final result for the two- loop box with two adjacent legs off-shell is: 

Jl m = J (jldwiTi-Wi)] {-s) d - 7 - w ^{-Mir\Miys{-ty« 

/(l, 1, 1, 1, w h w 2 , w 3 )f(l - w 2 , -w 3 , 1 - Wx, 1, w 4 , w 5 , w 6 ) (14) 

We have produced a representation for a two-loop integral by embedding representations for 
one-loop integrals into other representations. It is obvious that we can use the re-insertion 
method for writing representations for arbitrary multi-loop integrals. 

Representations for integrals with given kinematic scales can be used to derive, easily, 
results for simpler integrals where some of the scales are taken to zero. Let us, for illustration, 
derive the representation for the double-box with one-leg off-shell and the on-shell double- 
box. We must first take the limit M| — > in Eq. The term (M|) W5 is vanishing in 
this limit unless —>■ at the same time. If we use the Cauchy theorem, we find that 
the wanted limit is given by taking the residue of the integrand at w 5 = 0. Therefore, the 
double-box with one-leg off-shell is: 



-1 



(2m] 
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Jl m = t^-tt / [X[dw i T(-w t ))(-s) d - 7 - w ^{-Ml) w ^-tys 



\i=l 



/(l, 1, 1, 1, wi, w 2 , io 3 ).f(l - w 2 , -w 3 , 1 - wi, 1, iu 4 , 0, w 5 ) (15) 
Similarly, the on-shell double-box is given by: 

Jt m = ^4 / (jidwiTi-w^ (s) d -^(-tr* 

/(l, 1, 1, 1, wi, w 2 , w 3 )f(l - w 2 , -w 3 , 1 - Wi, 1, 0, 0, w A ). (16) 

Note that, in order to have a finite limit for a vanishing kinematic scale, the Mellin-Barnes 
representation should have terms in the integrand which behave as T(—g(wi ))M^ Wi \ These 
emerge naturally in non-trivial representations due to Eq. 01 
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III. ANALYTIC CONTINUATION 



In the previous Section we have noted that the Mellin-Barnes representation of a loop 
integral is valid if appropriate poles of the Gamma functions lay separated on the right 
and left of the integration contours. This condition guarantees the equivalence between the 
Mellin-Barnes integral and the original denominator in the loop integral. It often occurs that 
this condition can not be satisfied for values of the dimension parameter and the powers of 
the propagators which are relevant for a realistic application. Let us consider, as an example, 
the one-loop box representation of Eq. El in the case of powers of propagators set to unity 
and the dimension in d = 4 — 2e, 



'2m 



, dw 1 dw 2 dw 3 T(-w 1 )T(-w 2 )T 

{2m) 6 

xr (i + w 3 ) r (i + Wl23 ) r (-1 
-M?) W1 (-Miy 2 (-tr( 



-w 3 ) 
e - w 13 ) r 

\— 2— e— u>i23 



F(2 + e + w 123 ) 



T(-2e) 
-1 - e - w 23 ) 



(17) 



This representation can only be valid for values of e different from zero. For example, if 
we choose the contour C = {Re(wi) = — 0.1, Re(w 2 ) = — 0.2,Re(ws) = —0.3} we find that e 
should be in the interval —1.4 < e < —0.6. 

We can use the Cauchy theorem to obtain a representation in terms of a sum of contour 
integrals, which is valid at e = and the 1/e poles appear explicitly. The key point is that 
if the value of e is chosen outside the allowed region, for example e = 0, some of the poles 
of the Gamma functions will be on the wrong side of the contours (see Fig. Q). To recover 
the original result, we must compare the position of the poles relative to the contours for 
values of e inside the allowed region and the point outside. Depending whether the poles 
crossed the contours from left or right, we should correct the representation by adding or 
subtracting the residue of the integrand on those poles. This simple observation can be easily 
cast into general purpose algorithms which extract the poles in e of arbitrary Mellin-Barnes 
representations. 



Im(w) 



Im(w) 



e=0 



Re(w) 



Re(w) 



FIG. 1: On the left picture, all poles which originate from the same Gamma function are positioned 
either to the left or to the right of a contour of integration. On the second picture, we take e — > 0, 
and some of the poles cross to the other side of the contour. To recover a valid representation for 
the loop-integral, these poles sould be isolated using the Cauchy theorem. 

Let us choose a value of e inside the allowed region and observe the position of the 
poles with respect to the contours for some e-dependent Gamma functions in our example. 
The poles of T(— 1 — e — Wi 3 ) and T(— 1 — e — w 23 ), at — 1 — e — Re(w 3 ) + n > Re(u>x) 
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and — 1 — e — Re(u>3) + n > Re(w 2 ), are situated to the right of the wi and W2 contours, 
respectively, for all non-negative integers n. Let us now take e — > 0. We observe that the first 
poles, for n = 0, moved to new positions, — 1 — Ref/u^) < Re(wi) and — 1 — Re(iU3) < Re(w2), 
which are on the left side of the contours. The poles for n > remain on the right side. 
Therefore, the representation is not valid at e ~ 0, since the contours separate poles which 
originate from the same Gamma function. To obtain a valid expression for the integral at 
e = 0, we need to use the Cauchy theorem to isolate the crossing poles (n = 0). 

Smirnov, in Ref. |4l| . provides a number of pedagogical examples where this is done. The 
approach of Smirnov is to deform the multiple contours, using the Cauchy theorem, so that 
the "offending" poles never cross the contour of integration. Tausk, in Ref. [28j], describes 
in detail this task as well, using as an explicit example the Mellin-Barnes representation for 
the two-loop cross box integral. The technique of Tausk is different from the one of Smirnov; 
the main idea is to account for the poles that end up on the "wrong" side of the contour 
in a progressive way, as they cross the contours when changing continuously the value of 
e. In this paper, we have followed the approach of Tausk. We would like to recommend 
that the interested reader studies the pedagogical example of Ref. [Hf . Here we discuss our 
implementation of the algorithm which is depicted in the flow diagram of Figure El 

Consider the Mellin-Barnes representation 1Z of a loop integral. The integrand of TZ de- 
pends initially on Gamma functions. The application of the analytic continuation algorithm 
requires the iterative evaluation of residues. These may also give rise to Psi functions, which 
have poles in the same positions the Gamma functions do, and are treated identically for the 
purposes of analytic continuation. The contours C of the representation are straight lines 
parallel to the imaginary axis, and are chosen so that the real parts of the arguments of 
all Gamma functions are positive. The same condition determines an allowed region for e, 
which we will assume that includes the point e = eo- For the sake of clarity, we will consider 
the case eo < 0, but the method proceeds analogously when eo > 0. 

The first step is to find the maximum value e = ei, with ei < 0, for which the represen- 
tation holds. At this value, one or some of the positions of the poles run into a contour. If 
the set of these poles, V, is empty, i.e. no poles cross the contours when moving e from eo 
to 0, then the integral can be safely expanded around e = 0. So, let us focus on the more 
interesting case in which V is non empty. Here we have to distinguish between e\ = and 
ei < 0. 

If £i = 0, at least one pole lies on the contours for the terminating value of e. Then, the 
representation has an unregulated singularity, and it is not feasible to perform an e expansion 
and evaluate the contour integrals numerically. To remedy the situation we must change the 
contour of integration; we use the Cauchy theorem and displace one or more of the contours 
parallel to the imaginary axis. Then, we attempt afresh the analytic continuation; the poles 
should now hit the contour at e\ < 0. If this does not happen, due to an unfortunate choice 
of the contour we displaced, we simply iterate the procedure until we succeed. 

When €\ < 0, we have to determine the residues in the poles that cross the contour. If 
there is only one pole in V the situation is simple. Let us look back at the specific example 
of the one-loop box. By applying the previous steps of the algorithm, we find that the first 
pole to cross a contour is at — 1 — e — w 2 — = 0. We then choose one of the multiple 
integrals, e.g. the W2 integral, and subtract the residue of the integrand on this pole, at 
W2 = — 1 — e — W3, to the representation: 

K(e = e ) = TZ(e 1 ) - RaWI^-i-*-*, ( 18 ) 
The residue term in the above expression is a new multiple Mellin-Barnes integral with one 
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Solve the system of equations 
{U l =n l , Ui <EV}, 
for the MB variables: 
Wi = wf 




n = n + Y J] k 3 

eo = £1 



FIG. 2: Analytic continuation algorithm 
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integration variable less than in 1Z. In our example, the pole crosses the w 2 contour from 
the right side and we subtract its residue from the original representation. If a pole crosses 
the contour from the left, we should add its residue. At this point, we have achieved the 
analytic continuation in e from eo to e\. We now repeat the algorithm setting as initial value 
Eq = E\, and perform the analytic continuation to the next value of e where a new crossing 
occurs. We proceed iteratively, until we find all poles which cross the contours of integration 
in the terms of Eq. EH and its descendants. 

The algorithm is more complicated if many poles cross simultaneously the contour 
at e = €\. In our one-loop box example, if we had chosen a symmetric contour C = 
{Re(wi) = — 0.2, Re(w 2 ) = — 0.2,Re(u>3) = —0.2} then we would have two poles crossing at 
the same time: V = { — 1 — e — W13 = 0, — 1 — e — W23 = 0}. One way to deal with this situ- 
ation is to displace the contours in such a way that the degeneracy is removed, i.e. introduce 
asymmetric contours from the beginning. Another alternative is to take combined residues 
in all the poles. To correctly account for all the terms in this case, we can imagine that the 
contours are only slightly displaced from their current position. Then the poles will hit them 
one at the time and we can take residues successively. Finally, we can restore the original 
position of the contours without affecting the position of the poles. 

There is, however, a caveat in taking multiple residues. Let us consider the case of the 
one-loop box, where the poles Pi = — 1 — e — W13 and P 2 = — 1 — e — u> 23 hit the contour at the 
same time. We take residues in sequence, considering first P\. We make the "unfortunate" 
choice to take the residue on the W3 integration variable, and not on w\. The analytic 
continuation yields: 

K(e ~ e Q ) = K(e ~ ei ) - Rfis^)^.^ 

Then we consider the crossing of the P 2 pole on each of the two terms on the right hand 
side of the above expression. We find no problems in taking the residue of the first term. 
However, the position of the P2 pole for the second term is now written: 

P 2 = 1 - e - w 2 - w 3 L 3 =-i-e-™i = Wl - W2 ' 

The pole P 2 is now independent of e and gets nailed to the contour. This situation is similar 
to the one we encountered for e\ — 0, and we can deal with it in exactly the same way by 
displacing the contours. 

We apply the algorithm recursively, to all the terms which are generated by adding the 
residues of crossing poles to the original parameterization. At the end of this procedure, 
we obtain a sum of contour integrals plus some terms that do not contain any remaining 
integral due to taking residues on all integration variables. The method assures that the sum 
of all these contributions equals the original loop integral when e ~ 0. So we can expand in 
a power series of e. 

In this Section, we have presented the algorithm assuming that only the dimension pa- 
rameter e is required to regulate the Mellin-Barnes parameterization. However, the method 
proceeds identically for integrals where multiple analytic regulators are required. This could 
be the case in non-covariant gauges or in integrals with irreducible numerators [30] or addi- 
tional positive powers of propagators |29j ]. 

The actual implementation of the algorithm in a set of routines in both MAPLE and 
MATHEMATICA is very efficient, and allows to perform the analytic continuation in the 
regulator parameters in a few seconds, for all the Mellin-Barnes representations that we have 
studied in this paper. 
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IV. MELLIN BARNES REPRESENTATIONS FOR TENSOR INTEGRALS 



Mellin-Barnes representations, with the exception of the massless two-loop diagonal box 
topology have been traditionally employed for the evaluation of master integrals. In 
this Section, we introduce an efficient generalization of the method to tensor integrals; our 
method does not require any reduction techniques to master integrals. Common reduction 
methods produce extremely large algebraic expressions which may also have spurious sin- 
gular denominators. Here we derive efficient representations for one-loop tensor integrals, 
which are amenable to numerical integration. We can derive analogous representations for 
multi-loop integrals by using the re-insertion method of Section |Hj 

We consider a generic one-loop tensor integral of rank m with n external legs: 

r d d k k^...k^ m 
n ' m ~ J i^l \{k + gi y - ml] \{k + qtf -ml}--- \{k + q n f - m% (i9) 

where, 

9i = 0, <&• = £>, (20) 

i=l 

and pi are external incoming momenta. We first introduce Feynman parameters for the 
denominator. We obtain, 

I -T(n) [** f(f\dx)s(l-Tx) (21) 

«7T2 J \ i=1 / V i=l — 



where 



\(k + py- Af 

E^> (22) 



i=i 

and 

n n 

\2 



A = X i m i + EE X i X J [- (li ~ 1j) \ ( 23 ) 
i=l j=2 i<j 

Now we perform the usual shift in the loop momentum, k = K — P, and we obtain, 

= !■(„) /(n^) (A ,_ } A) „ ■ (24) 

where, we denote 

|^rpm-r jtW'-'^ = ^ . . . K Hr pN r+l . . . p« m ; (25) 

{jl,-,jm} 

{ji, ■ ■ -,jm} e permutations(l, . . . ,m). 

We find the standard Feynman representation for the generic tensor one-loop n-point func- 
tion, by integrating the loop momentum: 



X 



/(n^^i-E^^i 

E r ( W ~f ~S) {AP m - r } [W AS, (26) 



22 
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where A r = for r odd, and A r = g^ 1 ^ 2 ■ ■ ■ gVr-wr] f or r eV en. We observe that the sum, 
in the second line of the equation, is a polynomial in the Feynman parameters, with tensor 
coefficients: 

E T(jl ~J~^ {ArP-f 1 '-'^ Ai = £ C7(W; W})^" 1 • • -C- 1 (27) 

r<m Z {^i>0} 

For the scalar integral with unit powers of propagators, this sum is T(n — |). 

We will use Eq. |2H1 to derive Mellin-Barnes representations for the tensor loop-integral. 
We would like to exploit the fact that the Feynman representation of the scalar and tensor 
integrals have the same denominator: 

1 

In Section |Hj we have seen that we can obtain a Mellin-Barnes representation for a loop 
integral by decomposing the denominator of the Feynman parameterization with Eq. |3J We 
should anticipate that the Mellin-Barnes representations of tensors and scalars are, therefore, 
very similar. We can verify this intuition, if we consider one generic term in the polynomial 
of Eq. |23 and follow the steps of Section |TT1 Decomposing the denominator with Eq. El we 
obtain 

1 . . . rfVn — 1 

Yldw.ri-w^-s^i-s^A T(n - ~ + wi2...( a _i))a^ ^ • • -^28) 

where a is the number of kinematic scales in the integral. The factors xf 1 ^ are universal 
for all terms in Eq. (23 and originate from the dependence of A on the Feynman parameters. 
We now proceed to integrate out the Feynman parameters using Eq. [7| This integration 
yields Gamma functions 

x ^ )+n - 1 r(A + vi) = r(i + A)(i + A)(2 + A) • • • {n - 1 + A) (29) 

Rewriting all Gamma functions as above in Eq. 123 we obtain, for the tensor integral, the 
Gamma functions of the scalar integral multiplied by simple factors with Mellin-Barnes 
integration variables. 

The result for the one-loop tensor integrals depends on the topology and the rank of the 
tensors in a rather complicated manner. However, we find that the representation is of the 
following form 

In,m = ^ja-l" / (n^r(-^)(- S 0^ 

Xh^( Wl ,...,W a ^), (30) 

The first line of the above representation contains all the terms which already appear in 
the analogous representation of the scalar integral with unit powers of propagators. The 
function is a polynomial in the Mellin-Barnes variables, with tensor coefficients in terms 
of the external momenta. This is very important for our strategy to evaluate tensor integrals 
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and Feynman diagrams. The polynomial in the numerator is smooth, and it does not affect 
the analytic continuation in e. We can therefore perform the continuation collectively for all 
tensor integrals and diagrams of the same topology 

We organize the evaluation of tensor integrals and diagrams of a topology as follows. 
We first find the Mellin-Barnes parameterization of the scalar integral of the topology, and 
multiply the integrand with a template function h({wi}) which we assume depends on all 
Mellin-Barnes variables. Then, we apply the analytic continuation algorithm of Section 11111 
keeping h general. We expand in e, and create numerical programs for the evaluation of the 
expansion for any smooth function h. This part of the evaluation needs to be performed 
only once for a given topology. Then we must identify the polynomial for the evaluation of 
the specific tensors or diagrams that we are interested in. We have written FORM |43[ and 
MATHEMATICA programs which perform the steps that we described in this Section, and 
derive the explicit functional form of the polynomials in Eq. EH We then write numerical 
routines for their evaluation, and link them to the general routines for the e expansion of 
the integrals of the topology. 

This approach is efficient for the application of our technique to tensor integrals. It turns 
out, however, that the expressions for the polynomials are quite long for high rank tensors. 
We have observed that we can reduce the size of the expressions substantially with a simple 
modification. The terms A a in Eq. EHlare lengthy, and when we integrate out the Feynman 
parameters they give rise to large expressions in 

fcM of Eq. EH We rewrite Eq. M 

m a 

different way, 



:-«•/ (ft**) 



x e r( ":r i:) ' {A.r-f' "- 1 (3D 

r<m A2 A 2 

Then, we introduce a Mellin-Barnes decomposition for each denominator: 

1 



of 
Me 



and integrate out the Feynman parameters. At the end of this procedure, we obtain a sum 
Y + 1 Mellin-Barnes integrals for a tensor of rank m. All integrals correspond to the 
lin-Barnes representation of the scalar integral with shifted dimension d — > d + r. Each 
of them has a different polynomial in the numerator. However, these are substantially 
shorter than the polynomial in Eq. EDI 

It is worth to note that, with our method, we never introduce spurious singularities. The 
evaluation of the higher rank tensors and diagrams in gauge theories is very efficient, since 
we never require the manipulation of lengthy expressions. The polynomials can be 
large, however, they require a minimum amount of handling before generating numerical 
codes for their evaluation. 

Here we have presented the derivation of efficient representations for one-loop integrals. 
The method may be applied to multi-loop integrals using the one-loop tensor integral repre- 
sentations as building blocks. The analogous functions of in multi-loop tensors contain 
Gamma functions in the denominator. However, they are smooth and maintain the prop- 
erties which are required for a collective evaluation of all tensors of a topology. Explicit 
applications with multi-loop tensor integrals will be examined in future work. 
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V. NUMERICAL EVALUATION OF MELLIN-BARNES REPRESENTATIONS 



After analytic continuation, the contour integrals can be safely expanded in power series 
of e. The coefficients of the expansion are also contour integrals which, in general, contain 
in their integrand Gamma functions, Gamma function derivatives, simple logarithms, and 
powers of the kinematic parameters. In earlier works, these integrals were computed an- 
alytically. The standard approach is to evaluate them by means of the residue theorem, 
i.e. summing an infinite number of residues and expressing the power series as logarithms, 
polylogarithms and harmonic polylogarithms. 

This method is cumbersome for the evaluation of integrals with more than one Mellin- 
Barnes variable. In all cases presented in the literature, integrals with more than one 
dimension are reduced to one-dimensional integrals with a variety of clever, however not 
general, tricks. Typical procedures include direct integration in simple situations, the appli- 
cation of Barnes lemmas (27, 2^, 41 1, and the expansion in Laurent series around an integer 
value of one of the variables. In a few careful inspection of the multi-dimensional 

integrals reveals exact cancellations between different components which yield significant 
simplifications. 

Many integrals which are relevant for physical applications, such as the massless double 
box with two external legs off-shell and one-loop tensor integrals with six external legs, 
after the analytic continuation in e, are expressed in terms of hundreds of Mellin-Barnes 
components; some of them have a very high dimensionality. The analytic methods are 
not suitable for such problems. In our method, we choose to evaluate the MB integrals 
numerically by direct integration over the contours. 

The procedure we follow is, conceptually, straightforward. In practice, the implementa- 
tion of the numerical evaluation requires a significant work in automating the book-keeping 
of the various terms. Tensor integrals require significantly more complex book-keeping due 
to large expressions that multiply the Gamma functions in the integrand of their Mellin- 
Barnes representation. Our computer programs are organized in the following way. We first 
use MAPLE, FORM Q, and MATHEMATICA routines to derive Mellin-Barnes represen- 
tations for the integrals that we need to compute. We then perform the analytic continuation 
of Section 11111 and we expand the integrands in powers of e using MAPLE and MATHE- 
MATICA. As the next step, we translate the integrands into FORTRAN routines which 
evaluate them as complex quantities. 

There are many options to calculate numerically the sum of the Mellin-Barnes integrals 
which emerge from the routines for the e expansion. For example, we can combine all con- 
tributions in a single common integrand, so that large numerical cancellations take place 
before integration. However, the finding of the peaks and the adaptation of the integra- 
tion routines become less efficient. Following a diametric approach, we could integrate all 
components separately and compute their sum at the end. Then, the adaptation of the 
numerical integration is optimal, however, the numerical precision is sensitive to cancella- 
tions. In practice, we follow a hybrid approach and group together integrals with equivalent 
contours. We have achieved a reasonable efficiency and speed for all the computations that 
we present in this paper. Since we aim to keep our integration method general, we do not 
attempt any special simplifications by, for example, applying Barnes lemmas or exploiting 
symmetries of particular contours. 

For the numerical integration, we maintain the contours as straight lines parallel to the 
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imaginary axis. We map them onto a real interval with a simple transformation, 
1 r c+ioc 1 r 1 d\ ( ( \ W 

aS L. *" W = 5F I W^> f { c ~ 1 ln (i^a) ) (32) 

In all cases that we have considered, the integrands vanish rapidly when the integration 
variables take values away from the real line (A ~ |). The Mellin-Barnes integrals do 
not present problematic numerical instabilities, and converge rather fast. We perform the 
multidimensional integrations using the Cuhre and Vegas routines of the CUBA library j4^| ■ 
To gain some speed, we integrate the real and imaginary parts of the integrals concurrently, 
using the same grid. The grid is adapted to the peaks of the real part, however, the quality 
of the integration for the imaginary part is only mildly affected. On the other hand, there 
is another subtle point concerning the convergence of the imaginary parts. As they always 
start an order later in the expansion in e, the analytic expressions for them would be simpler 
than the corresponding real ones. When integrating numerically this relative simplicity is 
reflected by a faster convergence when compared to the real pieces at the same order. 

The numerical results depend on the values for the kinematic scales of the integrals. 
With our approach, we use the same expressions for the evaluation of the integrals in all 
kinematic regions. As we explained earlier, the only analytic continuations in kinematic 
variables that we must perform, is for simple functions of the type (— Sy)* and log(— sy). 
These are implemented trivially, in a branched format, in numerical programs. In this 
paper, we present multi-loop integrals for fixed kinematic parameters. However, we have 
arranged our routines so that a numerical integration over the phase-space can be performed 
simultaneously with the Mellin-Barnes integrations. 

A problematic case for our numerical algorithms is the evaluation of loop integrals with 
massive internal propagators in non-Euclidean kinematic regions. Typically, Mellin-Barnes 
represenations for such integrals exhibit a very slow damping at ±zoo in physical regions. 
The mapping of the integration region to a finite interval that we have considered in Eq. E21 
is not adequate, and contour deformations together with more sophisticated algorithms for 
the evaluation of oscillatory integrals are required. We defer the study of such situations for 
a future project. 

We have achieved a full automatization of the code generation for the numerical evalua- 
tion. This has allowed us to apply our method to diverse problems. Our approach can be 
improved and refined further in the future by considering, for example, a clever analysis of 
the terms before integration in order to avoid numerical cancellations or more sophisticated 
mappings which can smoothen the peak structure of the integrands. However, as we will 
see in the rest of the paper, our first, naive, implementation is already sufficient to obtain 
accurate results for complicated integrals in a reasonable time. 

In what follows we present results for one, two, and three loop integrals. Some of the 
integrals are known analytically, and we use them to verify our numerical programs. We also 
present new loop integrals, which would be extremely tedious to calculate with traditional 
methods. 

VI. RESULTS 

In this Section we present results for one, two, and three-loop integrals that we obtain 
with our method. For each integral, we present the coefficients of its expansion in e = 2 — ~ 
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through the finite term: 



J 



J \ \ mi J rij propagator^- 



r 



(33) 



where r is the depth of the leading pole in e. We estimate the errors associated with our 
results for q, by adding, in quadrature, the integration errors of all contributions from the 
residue decomposition. 

A. One loop hexagon: the scalar integral 

The first integral we consider is the massless one-loop hexagon, depicted in Figure 0J 
For massless external legs, we derive an eight-dimensional Mellin-Barnes representation for 



this integral. Due to the high dimensionality, the analytic continuation in e involves the 
crossing of the contours of integration from many poles, and the Mellin-Barnes expression 
after the continuation consists of about two hundred terms. This example demonstrates the 
advantage of the automatic algorithm we described in Section IIII1 since the book-keeping is 
done automatically and our routines perform the e expansion in fractions of a minute. 

After the analytic continuation and the expansion in series of e, only integrals with low 
dimensionality contribute to the poles and finite pieces of the series. The integration over 
Feynman parameters produces a denominator factor l/r(— 2e), which contributes factors 
of e in the series expansion. The eight-dimensional component, in which no residues have 
been taken, drops out from the finite part. In addition, most of the residues in the analytic 
continuation do not give rise to poles in e; a few components only develop a singularity that 
compensates the e in the overall factor. For the scalar hexagon, we find integrals with up to 
three dimensions which contribute to the poles and the finite term. This mechanism, which 
reduces the dimensionality of the representation in the leading coefficients of the e expansion, 
is not specific to the hexagon topology, but is typical of Mellin-Barnes representations of 
scalar multi-leg integrals. 

In Table ITD we present explicit results for the scalar hexagon in phase space points cor- 
responding to the physical region for the 2^4 process. To check the efficiency of our 





FIG. 3: The hexagon topology 
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integration routines in all the phase-space, we have sampled hundreds of points. Here, we 
only show representative results in the phase-space points of Table |U These points were 
generated starting from explicit values for the momenta in four dimensions, so they have a 
vanishing Gram determinant. 



point 


S12 


•S23 


S34 


•S45 


S56 


si6 


S123 


S234 


S345 


Pi 


1 


-0.232033 


0.096793 


0.025066 


0.465569 


-0.015783 


0.63356 


-0.219996 


0.336498 


P2 


1 


-0.056339 


0.054104 


0.111985 


0.583564 


-0.038557 


0.74554 


-0.191061 


0.260063 


P 3 


1 


-0.336912 


0.099306 


0.228826 


0.333228 


-0.086216 


0.64116 


-0.286003 


0.580452 


Pa 


1 


-0.310819 


0.012151 


0.020687 


0.561466 


-0.349136 


0.58216 


-0.318051 


0.114844 


P5 


1 


-0.07491 


0.048452 


0.279727 


0.314164 


-0.426813 


0.73852 


-0.449749 


0.532752 


P& 


1 


-0.447378 


0.10597 


0.052091 


0.277082 


-0.484988 


0.32942 


-0.400468 


0.170418 


Pi 


1 


-0.081759 


0.143257 


0.139358 


0.183411 


-0.544412 


0.81512 


-0.361942 


0.31596 


P 8 


1 


-0.115263 


0.054869 


0.154172 


0.209043 


-0.569659 


0.63722 


-0.477605 


0.410718 


P9 


1 


-0.331123 


0.011523 


0.033009 


0.252287 


-0.635452 


0.38202 


-0.441643 


0.309497 


Pw 


1 


-0.297484 


0.010637 


0.05663 


0.345614 


-0.88129 


0.68562 


-0.541753 


0.083132 



TABLE I: Explicit values of the invariants for some of the phase space points used for the evaluation 
of the hexagon topology. The points correspond to the physical region of the 2 — > 4 scattering of 
massless particles. 

As can be seen from Table |Hj our method is perfectly capable of evaluating multi-scale 
integrals, such as the scalar hexagon, reliably in all phase-space. The results in the Table 
correspond to a fixed number of evaluations, thus the difference in the relative errors for 
different points. It took a couple of minutes per point on a 2.8GHZ CPU to evaluate them. 

B. One-loop hexagon: tensor integrals 

The most severe problem in evaluating multi-leg amplitudes, is the appearance of tensor 
numerators in the loop integrals. Our method, provides a powerful tool for their evaluation. 
The scalar hexagon, which we discussed before, provided a nice test-ground for the analytic 
continuation algorithm and the numerical integration strategy. Here, we consider tensor 
integrals for the same topology, in order to demonstrate the potential of the method for 
computing realistic quantities that arise in gauge theories. 

Using the procedure described in Section HVl we perform the analytic continuation for an 
arbitrary tensor. At variance with the scalar case, now we find components with more than 
three MB variables contributing to the finite pieces in e. For instance, for a typical rank 
three tensor, integrals with all eight Mellin-Barnes variables are required. This calculation 
provides stringent tests on our implementation, due to the complex book-keeping of the 
tensorial terms in the Mellin-Barnes representations, and the efficiency of our routines to 
evaluate the integrals with high dimensionality. 

In Table 11111 we show results for tensors up to rank 6 evaluated in a symmetric point 
{sij = Sijk = —1} in the Euclidean region. For these examples, we contracted all the loop 
momenta in the numerator with the same external momentum. All such tensor contractions, 
with odd rank, should vanish; we reproduce this result. On the table, we also show results 
obtained by reducing the tensors using the program AIR We were not able to reduce 
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Point 


C2 


Cl 


co 


Pi 


- 24550.802400 (0) 
+ %■ 0.000000 (0) 


- 129377.216040 (5) 

- i • 69329.856499 (5) 


249435. (40.) 
- »• 365517.9 (3) 


P2 


- 10666.332700 (0) 
+ %■ 0.000000 (0) 


42307.989040 (3) 
- i • 18370.333000 (3) 


34595. (66.) 
- * - 53139.1 (8) 


P3 


- 1153.682480 (0) 
+ i- 0.000000 (0) 


3091.953473 (3) 
- i ■ 2868.081530 (4) 


2292. (2.) 
- i- 6531.15 (6) 


Pa 


- 48753.897600 (0) 
+ i- 0.000000 (0) 


256890.405260 (4) 
- i ■ 152982.841000 (4) 


- 432690. (5.) 

- *■ 806442.61 (4) 


P5 


- 2502.711680 (0) 
+ i- 0.000000 (0) 


- 10214.467200 (2) 

- i • 7425.550140 (3) 


8934. (16.) 
- %■ 30275.25 (6) 


P(> 


- 3857.953670 (0) 
+ i- 0.000000 (0) 


11372.775324 (2) 
- i • 12047.807600 (3) 


6715.7 (9) 
- i- 35512.41 (2) 


Pi 


- 2078.190700 (0) 
+ i- 0.000000 (0) 


6541.063310 (1) 
- i ■ 6138.816950 (2) 


1413. (2.) 
- i • 18965.48 (6) 


Ps 


- 3356.184660 (0) 
+ i- 0.000000 (0) 


13032.698020 (3) 
- i • 10343.407301 (2) 


8558. (6.) 
- i- 40209.00 (3) 


P9 


- 45864.213500 (0) 
+ %■ 0.000000 (0) 


281098.336280 (4) 
- i ■ 144019.061999 (3) 


649614. (5.) 
- %■ 882892.34 (2) 


P10 


- 21516.252000 (0) 
+ %• 0.000000 (0) 


97739.299610 (3) 
- i ■ 67551.061100 (3) 


99008. (13.) 
- %■ 306919.34 (2) 



TABLE II: The scalar massless hexagon evaluated in some physical points of the phase space. 
The errors of the numerical integration are quoted in parenthesis, and they affect the last figure/s 
respectively. 



the high rank tensors to master integrals in a closed analytic form. The reduction was only 
possible by substituting the kinematic invariants to their numerical value at this particular 
phase-space point before solving the system of IBP equations. We have also made additional 
comparisons with the reduction method, for other tensor contractions and a different point 
inside the physical region. 

As in the scalar case, we sampled several phase space points in the physical region, 
calculating a variety of contracted tensors at each rank. In Table IIVI we list the cases we 
considered, and in Table we show explicit results for the contraction q2 ■ k q2 ■ k q% ■ k ■ 
kq$ ■ kq e ■ k in the same phase space points that we considered for the scalar hexagon. 
We find that our routines perform extremely well, including the highest rank cases that we 
considered here. 

We should note that, due to the way that we organize the evaluation of the tensors in 
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Tensor 


MB 


AIR 


<?2 • k 


(0.000000 (0) + i-0.000000 (0))e- 2 + 
(0.00000 (1) +Z-0.00000 (l))e _1 + 
(0.000 (2) +i-0.00000 (1)) 





q 2 -kq 2 -k 


(0.500000 (0) +i-0.000000 (0))e- 2 + 
(0.71139 (1) +Z-0.00000 (l))e _1 + 
(0.0947 (8) +i-0.00000 (1)) 


0.5 er 2 + 
0.711392 e _1 + 
0.094845 


q 2 - kq 2 ■ kq 2 - k 


(0.000000 (0) +i-0.000000 (0))e~ 2 + 
(0.0000 (1) +Z-0.00000 (l))e _1 + 
(0.00 (2) +i-0.02 (2)) 





q 2 ■ kq 2 ■ kq 2 ■ kq 2 ■ k 


(0.125000 (0) +i-0.000000 (0))e- 2 + 
(0.3861 (1) +i-0.00000 (l))e~ x + 
(0.65 (3) +i-0.02 (3)) 


0.125 e~ 2 - 
0.386181 e -1 + 
0.65996 


q 2 ■ kq 2 ■ kq 2 ■ kq 2 ■ kq 2 ■ k 


(0.000000 (0) +i-0.000000 (0))e- 2 + 
(0.00000 (9) +i-0.00000 (l))e- 1 + 
(0.04 (4) -i-0.02 (4)) 





q 2 ■ kq 2 ■ kq 2 ■ kq 2 ■ kq 2 ■ kq 2 ■ k 


(0.031250 (0) +Z-0.000000 (0))e- 2 + 
(0.12466 (6) +i-0.00000 (l))e- 1 + 
(0.27 (1) -i-0.00 (1)) 


0.03125 e" 2 - 
0.1246703 e~ l + 
0.279587 



TABLE III: Contracted tensors evaluated in a symmetric point in the Euclidean region. We 
compare with results obtained with the program AIR. 

Section II VI we do not anticipate that the evaluation of Feynman graphs, such as the ones 
in the one-loop six photon amplitude, are significantly harder than the six rank tensors that 
we presented. In addition, our routines allow a combined integration over the Mellin-Barnes 
variables and the kinematic invariants. 

C. On-shell planar double box 

Pi P4 



P2 P3 

FIG. 4: The massless double box. 
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Rank 


Contraction 


1 


92 • k 


2 


q 2 ■ kq 4 - k 
q 2 -kq 5 -k 


3 


rtn • Jc n<~> • Jc a a - Jc 

q 2 ■ kq^- kq 5 ■ k 


4 


q 2 ■ kq 2 ■ kq 3 ■ kq A ■ k 

y2 </4 ys "'Ho 
q 2 ■ kq 3 ■ kq 5 ■ kq 6 ■ k 


R 


q 2 ■ kq 3 ■ kq A - kq 5 ■ kq 6 ■ k 

t/2 y2 y4 96 ™ hk> 
k ■ k q 2 ■ k q 3 ■ k qi ■ k 


6 


q 2 -kq 2 -kq 3 -kq4-kq 5 -kq 6 -k 
k ■ k q 2 ■ k q 3 ■ k q^ ■ k q$ ■ k 

If . If If . If 
iXi fx, i\j rv hj rv 



TABLE IV: Tensor contractions evaluated for the hexagon topology 



We now compute two-loop integrals, and we consider first the massless planar double 
box with on-shell legs shown in Fig. 0] This integral is known analytically in Ref. j2?| . 
As discussed in Section |HJ it can be expressed as a MB contour integral in four complex 
dimensions. Performing the analytic continuation in e, we get terms where residues in all 
the variables have been taken (i.e. terms without any contour integral left); these contain 
poles through 1/e 4 . The terms through O(e ) require Mellin-Barnes integrals with one and 
two integration variables. 

In Figure El we show our results for the finite pieces of the on-shell double box in the 
physical region for the pi +p 2 — > p 3 +p^ process depicted in FigurelH We set s = (p± +p 2 ) 2 = 
1, and plot, as functions of t = (p 2 — p 3 ) 2 , the finite term Co and the comparison of our 
numerical calculation to the analytic result of |27| . 

As can be seen from the figure, it is possible to achieve an accuracy of a few parts in ten 
thousand, in a few seconds per phase-space point. The bigger error band for the imaginary 
part, reflects our naive strategy to integrate both real and imaginary part with the same 
grid. The coefficients for the simple and double poles in e, not shown in the figure, agree 
with the analytic calculation with even better accuracy, whereas the cubic and quartic poles 
are analytic expressions also in our case (as they do not involve any MB integration) and 
coincide with the ones in I27II. 



D. On-shell non-planar double-box 

Now we compute the non-planar two-loop box with massless external legs depicted in 
Fig. |U1 This integral has been computed analytically in Ref. |28| using a four-fold Mellin- 
Barnes representation. Instead of deriving a representation with the aid of the re-insertion 
method, we use the the representation in Ref. |28( for our numerical evaluation. This selec- 
tion allows us to check our algorithms for the expansion in e with the detailed decription 
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Point 


C2 


c\ 


Co 


Pi 


- 402.790580 (0) 
+ *• 0.000000 (0) 


1919.50 (1) 
- % ■ 1248.33 (2) 


3019. (55.) 
- i-6068. (50.) 


P2 


- 253.831615 (0) 
+ i- 0.000000 (0) 


876.7 (2) 
- %■ 676.4 (1) 


538. (56.) 
- i-2440. (50.) 


P3 


- 47.453378 (0) 
+ %■ 0.000000 (0) 


149.099 (6) 
- i ■ 145.128 (4) 


79. (2.) 
- i-485. (2.) 


Pa 


- 3876.804240 (0) 
+ %■ 0.000000 (0) 


23083.200 (2) 
- % • 12161.074 (2) 


53835. (6.) 
- % ■ 72507. (5.) 


P5 


- 52.452711 (0) 
+ i- 0.000000 (0) 


214.000 (2) 
- % ■ 158.220 (2) 


214. (2.) 
- z-663. (1.) 


P(> 


- 224.589492 (0) 
+ i- 0.000000 (0) 


887.7120 (5) 
- i- 701.8064 (5) 


888. (1.) 
- *■ 2796.9 (9) 


Pi 


- 54.532983 (0) 
+ i- 0.000000 (0) 


172.177 (4) 
- i ■ 161.007 (7) 


26. (2.) 
- z - 531. (2.) 


Ps 


- 47.254266 (0) 
+ i- 0.000000 (0) 


179.409 (1) 
- i ■ 143.862 (2) 


127. (1.) 
- z-558. (1.) 


P9 


- 865.507076 (0) 
+ i- 0.000000 (0) 


5612.7500 (2) 
- i- 2716.2898 (1) 


15300. (2.) 
- i ■ 17634. (1.) 


P10 


- 1318.087970 (0) 
+ i- 0.000000 (0) 


6160.7460 (9) 
- %■ 4137.4238 (2) 


7611. (2.) 
- i ■ 19356. (2.) 



TABLE V: Results for the rank six tensor qi ■ k q2 ■ k q% ■ k q^- k q$ ■ k ■ k for the hexagon topology 
at one loop. The phase space points are detailed in Table [I] 



in Ref. 28]. Furthermore, this representation does not impose the physical constraint be- 
tween the three invariants, s + t + u = 0, leading to a more complicated analytic structure. 
This is an additional test on our programs for evaluating integrals in phase space regions 
which are separated by complicated thresholds. The integral is, also, 1/e 4 divergent. Up 
to three-dimensional integrals contribute to the e expansion through O(e ). In Ref. [28], 
it was shown that, with clever manipulations, it is possible to reduce all multiple integrals 
to integrals with only one dimension. Keeping our routines general, we have chosen not to 
reduce the dimensionality of the multiple integrals and evaluate them numerically. 

In Figures [7| and |H1 we show our results for the finite piece of the non-planar box in the 
regions (i) u, t < and s = —t — u and (ii) u, s < and t = —s — u respectively, where 
s — (pi + P2) 2 , t = (j>2 — P3) 2 and s + t + u = 0. In case (i) we fixed s = 1 as reference and 
plot the results as a function of t, whereas in (ii) t = 1 is fixed and we change s within its 
allowed range. We show the comparison of our results with the analytic calculation of |28| . 
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FIG. 5: Results for the finite part of the planar double box in the physical region for a 2 — > 2 
process. On the two left panels we plot the real and imaginary parts (upper and lower plots 
respectively) of the finite term as a function of t for fixed value of s = 1. The estimated error of 
the numerical integration lies within the size of the points. On the right panel we show the ratios 
of the numerical calculation to the analytic results of f° r the same kinematics, the bands in 
this case are given by the error in the numerical integrations. 



In this case, the relative errors obtained by numerical integration are bigger than in the 
the on-shell double box. This is particularly noticeable in the imaginary part in both regions 
(i) and (ii), for which the relative error reaches 1-2%. However, this only happens for non- 
significant points of the phase-space in which the magnitude of the imaginary part is almost 
zero (notice that in both regions, the imaginary part of the finite term changes sign twice). 

The non-planar double-box has a complicated analytic structure in terms of the kine- 
matic invariants, since there is no Euclidean kinematic region for this integral. The two 
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FIG. 6: The cross box. 

kinematic regions that we study here, require complicated analytic continuations with tra- 
ditional methods. With our method, it is particularly simple to compute loop integrals in 
different regions of phase space. Our results demonstrate that the numerical integration 
of the contour integrals allows for a trivial analytic continuation in the invariants, giving 
correct and accurate results in all phase-space. 



E. Planar double-box with one leg off-shell 

Our next example is the planar double-box with one external mass. The first results 
for this integral were obtained in |22| using sector decomposition. Soon afterwards it was 
computed analytically in j3l|, and with differential equation method in (37l |. 

We derive a 5-fold MB representation for this integral in Section II. The integral is 1 /e 4 
divergent; after the expansion in e, we find that the double and simple poles and the finite 
term, contain three-dimensional Mellin-Barnes integrals. 

In Figure ITU1 we show the results for this integral in the region corresponding to the decay 
of a heavy particle p 4 — > P1+P2+P3 (Figure^), where we have fixed the mass of the particle 
v\ — s i23 = 1 and one of the invariants S13 = (pi + p 3 ) 2 = 3/10, and we consider values of 
S23 = {p2 +P3) 2 within the allowed range. Again, the relative errors are of the order of the 
few per mille, with highest values in the points in which the absolute value of the integral 
is small (in this case both the real and the imaginary components change sign). 

For this topology we also considered the kinematical region corresponding to a p\ + P2 — > 
P3 + P4 process being p 4 the massive momentum (see Fig. The corresponding results 

are shown in FigurelTTlfor fixed values of s = (p\ + P2) 2 = 1 and pi = 1/10 as a function of 
t = (P2 ~ Pz) 2 - Again we find very good agreement with the analytic results in both phase 
space regions. 



F. Planar double box with two adjacent massive legs 

We now make one further step for the planar double box and consider the case of two 
external adjacent masses. This integral has been evaluated for the first time, in some points 
in the Euclidean region, in j23|. 

As shown in Section II, it is possible to get a MB representation for the double box with 
two adjacent massive legs with six MB parameters. Again, after the analytic continuation, 
the effective dimensionality is reduced by factors containing gamma functions in the denom- 
inator. In the present example, only integrals in three dimensions are needed to compute the 
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FIG. 7: Results for the finite piece of the non planar box with seven propagators in the region 
u, t < and s = —t — u for s = 1 as a function of t. The left panels show our numerical 
results for the real and imaginary part (upper and lower plot respectively) and the right ones the 
corresponding comparisons with the analytic calculation of 



finite pieces in the series around e = 0. Quartic and cubic poles are completely determined 
by the pieces with no integrals left. 

We found complete agreement with the results quoted in j^. We also compared our 
results for some values of the invariants in the Euclidean region with an independent calcu- 
lation obtained with our own sector decomposition code. In Table IVII we show the results 
obtained by the two methods for two phase space points in the mentioned region. The 
invariants are defined as s = {p\ + P2) 2 , t = (j>2 — P3) 2 , Mf = p|, Mf = pi (see Figure H21 for 
the definition of the momenta). We find perfect agreement within the integration errors. 

We also present our results in the kinematical region corresponding to the process p\ + 
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s s 

FIG. 8: Results for the finite piece of the non planar box with seven propagators in the region 
u, s < and t = —s — u for t = 1 as a function of s. The upper and lower left panels show 
our numerical results for the real and imaginary parts, respectively; and the two right ones the 
corresponding comparisons to the analytic calculation of 28]. 

Vi ~~ > Ps + Pi shown in Figure El relevant, for instance, for the process pp — > W + W~ . In 
Figure l*HH we plot the values of the coefficients, for this configuration, up to the finite terms 
as a function of t in the case of two equal masses, Mf = Mf = 1/10, where the energy 
scale is fixed by s — 1. Figure ITU corresponds to the same kinematical region but for s = 1, 
Mf = 1/20 and Mf = 1/2. The error bars again lie within the points in the plot and are, 
for the finite pieces, typically less that 1% after a run of a couple of minutes per point on a 
2.8GHZ CPU. 
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(a) 



(b) 



FIG. 9: The two kinematical configurations considered for the double box with one massive leg. 
(a) the decay process P4 — > p\ + p 2 + P3 and (b) the scattering p\ + p2 — > P3 + p±. The double line 
corresponds to the massive particle, p^. 



Point 


MB 


Sector Decomposition 


s = t = Ml=M% = -1 


-0.25 e~ 4 + 
0.288608 e- 3 + 
-2.22227 (1) e" 2 + 
-6.577 (6) 
-10.15 (4) 


-0.25 e- 4 + 
0.2886 (2) e-3+ 
-2.222 (2) e" 2 + 
-6.577 (6) (T l + 
-10.15 (2) 


s=-l 
t=-l/2 
Mf=-7/W 
M$=-2/5 


-0.5 e- 4 + 
0.463887 e- 3 + 

-4.47418 (1) e" 2 + 
-17.46 (1) e _1 + 
-37.21 (9) 


-0.5 e~ 4 + 
0.4639 (5) e~ 3 + 

-4.474 (4) e" 2 + 
-17.46 (2) e- x + 
-37.21 (9) 



TABLE VI: Results for the double box with two legs off-shell for two points in the Euclidean region, 
we show results obtained with both the Mellin-Barnes technique and wit the sector decomposition 
method. 

G. The on-shell massless triple box 

At the three loop level, we start by considering the on-shell triple box (Figure IT3j) . im- 
pressively computed in an analytic way in [32] using the MB technique. 

Using the re-insertion method described in Section II, we get a 7-fold MB representation. 
The analytic continuation leaves us with contributions with up to 5 contour integrals for the 
constant pieces in the e expansion, and poles up to e~ 6 . Although the numerical integration 
is certainly more involved in this case, due to the depth of the poles in e which generates 
very complicated expressions for the finite terms, it is relatively easy to achieve a precision 
of ~ 1%, as shown in Figure ITrJl As in the two loop case we show results for the physical 
region corresponding to s > 0, t < with s = Jj>i + P2) 2 and t = (p 2 — Pz) 2 - We compare 
our results with the analytic computation of [32]. We used the MATHEMATICA package 
HPL |4^] for the evaluation of the one-dimensional harmonic polylogarithms. 
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FIG. 10: Results for the finite term of the double box with one leg off-shell in the physical region 
for the decay of a heavy particle, p^ — > p\ + P2 + P3 (Figure On the upper and lower left panel 
we plot the real and imaginary parts as a function of the invariant S23 = (j>2 +P3) 2 for fixed value 



of pi 



S123 



1 and S13 = (pi + P3) 2 = 3/10. On the right panels we show the corresponding 



ratios of the numerical calculation to the analytic result of 
are given by the error in the numerical integrations. 



for the same kinematics, the bands 



H. The triple box with one external mass 

As we did for the double box topology, we consider now a further step in complication 
for the triple box, adding a mass in the external state. This is, as far as we know, the first 
evaluation of a three-loop box with three mass scales. 

We obtained a MB representation with 8 variables. As usual, the dimensionality of the 
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FIG. 11: Results for the finite part of the double box with one leg off-shell in the physical region 
for a 2 — > 2 process with the massive leg on the final state (Fig. |Hp). On the upper and lower left 
panels we plot the real and imaginary parts of cq as a function of the invariant t = (jp2 — P3) 2 for 
fixed value of s = (pi +P2) 2 = 1 and p\ = 1/10. On the right panels we show the corresponding 
ratios of the numerical calculation to the analytic result of for the same kinematics, the bands 
are given by the error in the numerical integrations. 



problem is reduced after the analytic continuation and we found that, up to order e°, only 
up to six-fold integrals contribute. 

We present results for this novel integral in the physical region of the decay of a heavy 
particle, p 4 — > pi +p 2 +P3 (Figure [T7J) as we did in the two loops case. In Figure [TH] we show 
the the corresponding coefficients of the expansion in e as a function of the invariant S23 for 
fixed values of the other two variables: p\ — \ and S13 = 3/10. With the notable exception 
of the left-most point in the real part of the constant term, the error bars are all contained 
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FIG. 12: The double box diagram with two adjacent massive legs. 

in the plotted points. They are typically of the order of 1% or better for the finite parts 
after 50 minutes, in average, per point on a 2.8GHZ CPU. 

VII. CONCLUSIONS AND OUTLOOK 

In this paper, we have introduced a new method for the numerical evaluation of arbi- 
trary loop integrals. We first derive Mellin-Barnes representations by either decomposing 
directly the Feynman representation, or by using one-loop representations as building blocks 
for representations of multi-loop integrals. For a given representation, we find values of the 
space-time dimension and the powers of the propagators which yield a well defined repre- 
sentation. Then we use automated programs to perform the analytic continuation to values 
of the parameters for which the integral may develop divergences. Finally, we expand in e 
and compute the integral coefficients of the expansion numerically. 

The method is based on the earlier work of Smirnov ji?} and Tausk [i^. The novel fea- 
tures in our publication are: (i) the automation of the procedure for the e expansion, (ii) the 
numerical evaluation of the coefficients of the expansion, avoiding cumbersome summations 
of multiple infinite series and the analytic continuation in the arguments of polylogarithms, 
(iii) and the efficient generalization of the method to tensor integrals. We believe that these 
new features render our method suitable for the practical evaluation of multi-loop amplitudes 
in gauge theories. 

We have performed a number of explicit calculations to verify our algorithms and demon- 
strate the power of our method. We first recalculated complicated loop integrals which are 
only recently known in the literature. Namely, we evaluated the planar and cross on-shell 
double boxes, the planar double-box with one leg off-shell, and the triple planar box with 
on-shell legs. In all cases, we have found an excellent agreement between our numerical 
results and the known analytic results. The analytic results require complicated analytic 
continuations of polylogarithms to non-Euclidean kinematic regions. With our method, we 
were able to compute the integrals in all kinematic regions effortlessly. 

In this paper, we have presented results for loop-integrals which were previously unknown. 
At two-loops, we have computed the double-box integral with two adjacent legs off-shell for 
independent external mass-parameters. This is one of the most complicated two-loop box 
integrals that enters the evaluation of two-loop amplitudes for heavy boson pair production 
at colliders. At three-loops, we have computed the planar triple-box with one off-shell 
leg; this integral emerges in e + e~ — > 3 jets, or the production of a single heavy boson in 
association with a jet at hadron colliders at NNNLO in QCD. To the best of our knowledge, 
it is the first time that a two-loop box with four kinematic scales or a three-loop box with 
three scales are ever computed in all physical regions. 
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FIG. 13: Results for the double box with two adjacent masses in the physical region of a 2 — > 2 
process where the two particles in the final state are massive. We plot the values of the coefficients 
as function of t for s = 1, M\ = Mf = 1/10. 



32 



-7 



16 -OA -0.2 

t 





■0.6 -0.4 



FIG. 14: Results for the double box with two adjacent masses in the physical region of a 2 — > 2 
process where the two particles in the final state are massive. We plot the values of the coefficients 
as function of t for s = 1, M\ = 1/20 and Mf = 1/2. 
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FIG. 15: The on-shell triple box diagram. 



Many processes with six external legs are particularly important at the LHC At present, 
there has been no NLO calculation of a cross-section for a hadron collider process with six 
external states, due to the lack of efficient methods for evaluating loop amplitudes. In this 
paper, we applied our method to the evaluation of tensors through rank six for the hexagon 
topology. The purpose of this application was two-fold; first, we wanted to demonstrate that 
we were able to extend the method to tensor integrals economically, and second, to setup 
programs that are efficient for the evaluation of multi-scale one-loop amplitudes (e.g. six 
parton QCD amplitudes). We have found that the evaluation of the tensor integrals is not 
significantly more complicated than the scalar integral. The programs that we developed 
can be immediately used for the evaluation of the six parton one-loop QCD amplitudes. 

Our technical goals for the current study were to build the necessary programs, implement 
an efficient book-keeping platform, and to examine the scaling of computing intensity in 
applications with diverse features. We anticipate that our method can be improved and 
generalized even further, when we consider issues that were only briefly examined in this 
first study. For example, in preliminary investigations, we have found that it is possible 
to exploit further the Cauchy theorem. If we specify a kinematic region, it turns out that 
some of the numerical integrals can be approximated very accurately by the sum of a finite 
number of residues. We expect that we can improve substantially our efficiency if we replace 
some of the integrations which involve kinematic scales with appropriate finite sums of such 
residues. In addition, we are investigating hybrid approaches, combining our method with 
ideas in Ref . . 

The applications of our method are numerous. We believe that it will be particularly 
useful for precision calculations of observables in collider physics. We plan to apply our 
techniques to key NLO and NNLO calculations for the LHC and to answer more formal 
questions on the perturbative behavior of gauge theories in the near future. 
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FIG. 16: Results for the finite part of the planar triple box in the physical region for a 2 — > 2 
process. On the upper and lower left panel we plot the real and imaginary parts, respectively, of 
Co as a function of t for fixed value of s = 1. On the ri ght panels we show the corresponding ratios 
of the numerical calculation to the analytic result of |32j for the same kinematics, the bands are 
given by the error in the numerical integrations. 
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FIG. 18: Results for the triple box with one massive leg in the physical region of the decay of a 
massive particle, p^, — > p\ + P2 + P3 (Figure IT7|) . We show the results for fixed values of pi = 1 and 
S 13 = 3/10 as a function of the remaining invariant S23. 
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FIG. 19: Results for the triple box with one leg off-shell in the physical region for a 2 — > 2 process 
with the massive leg on the final state. We plot the real and imaginary parts as a function of the 
invariant t = (p2 — P3) 2 for fixed value of s = (p\ + P2) 2 = 1 and p\ = 1/10. 
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